package geo

import (
	"errors"
	"fmt"
	"math"

	"gitee.com/ymofen/gobase"
	"github.com/fogleman/delaunay"
)

type gridQuadrangle struct {
	tri1       [3]int
	tri2       [3]int
	commonEdge [2]int
}
type Quadrilaterals struct {
	quadrs        []gridQuadrangle
	triangulation *delaunay.Triangulation
}

// 借鉴python interp_values = griddata((x_points, y_points), values, (x_grid, y_grid), method='linear')
func Griddata(x_points, y_points, values, x_grid, y_grid []float64) (rval []float64) {
	rval = make([]float64, 0)
	trion, _ := Triangulate(x_points, y_points)
	Quad := Quadrilaterals{triangulation: trion}
	Quad.constructingQuadrilaterals()
	// Quad.printSelf()
	// tris := Quad.triangulation
	for index := range x_grid {
		p := delaunay.Point{X: x_grid[index], Y: y_grid[index]}
		quadIndex, err := Quad.findTriByPoint(p)
		if err == nil {
			q := Quad.getQuad(Quad.quadrs[quadIndex].tri1, Quad.quadrs[quadIndex].tri2)
			if Quad.quadrs[quadIndex].tri2 == [3]int{0, 0, 0} {
				v := calculateTByTri(Quad.triangulation.Points[q[0]], Quad.triangulation.Points[q[1]], Quad.triangulation.Points[q[2]], p, values[q[0]], values[q[1]], values[q[2]])
				rval = append(rval, v)
			} else {
				v := (calculateT(Quad.triangulation.Points[q[0]], Quad.triangulation.Points[q[1]], Quad.triangulation.Points[q[2]], Quad.triangulation.Points[q[3]], p, values[q[0]], values[q[1]], values[q[2]], values[q[3]]))
				rval = append(rval, v)
			}
			// fmt.Println(x_grid[index], y_grid[index], ":", tris.Points[Quad.quadrs[triIndex].tri1[0]], tris.Points[Quad.quadrs[triIndex].tri1[1]], tris.Points[Quad.quadrs[triIndex].tri1[2]])
		} else {
			rval = append(rval, -1)
		}
	}

	// var lngs, lats []float64
	// points := trion.Points
	// tri := trion.Triangles
	// for ti := 0; ti < len(tri); ti++ {
	// 	lngs = append(lngs, points[tri[ti]].X)
	// 	lats = append(lats, points[tri[ti]].Y)
	// }
	// triTofile(lngs, lats)
	// pointTofile(points, "all")
	return
}

func (this *Quadrilaterals) findTriByPoint(point delaunay.Point) (int, error) {
	for i, one := range this.quadrs {
		if this.isPointInTriangle(point, one.tri1) {
			return i, nil
		}
	}
	return -1, errors.New("此点不在任何三角形内")
}

func (this *Quadrilaterals) appendQuadr(i, j, commonA, commonB int) {
	tris := this.triangulation.Triangles
	if j == -1 {
		this.quadrs = append(this.quadrs, gridQuadrangle{
			tri1: [3]int{tris[3*i+0], tris[3*i+1], tris[3*i+2]}})
		return
	}
	this.quadrs = append(this.quadrs, gridQuadrangle{
		tri1:       [3]int{tris[3*i+0], tris[3*i+1], tris[3*i+2]},
		tri2:       [3]int{tris[3*j+0], tris[3*j+1], tris[3*j+2]},
		commonEdge: [2]int{commonA, commonB}})
}

func (this *Quadrilaterals) constructingQuadrilaterals() {
	this.quadrs = []gridQuadrangle{}
	tris := this.triangulation.Triangles
	for i := 0; i < len(tris)/3; i++ {
		hadFindConvexQuad := false
		for _, one := range this.findTrisByLine(tris[3*i+0], tris[3*i+1]) {
			if !hadFindConvexQuad && one != i && this.isConvexQuad(i, one) {
				//fmt.Println(i, one)
				hadFindConvexQuad = true
				this.appendQuadr(i, one, tris[3*i+0], tris[3*i+1])
			}
		}
		if hadFindConvexQuad {
			continue
		}
		for _, one := range this.findTrisByLine(tris[3*i+1], tris[3*i+2]) {
			if !hadFindConvexQuad && one != i && this.isConvexQuad(i, one) {
				//fmt.Println(i, one)
				hadFindConvexQuad = true
				this.appendQuadr(i, one, tris[3*i+1], tris[3*i+2])
			}
		}
		if hadFindConvexQuad {
			continue
		}
		for _, one := range this.findTrisByLine(tris[3*i+2], tris[3*i+0]) {
			if !hadFindConvexQuad && one != i && this.isConvexQuad(i, one) {
				//fmt.Println(i, one)
				hadFindConvexQuad = true
				this.appendQuadr(i, one, tris[3*i+2], tris[3*i+0])
			}
		}
		if !hadFindConvexQuad {
			// fmt.Println(i, tris[3*i+0], tris[3*i+1], tris[3*i+2], "未找到凸四边形")
			this.appendQuadr(i, -1, 0, 0)
		}
		// fmt.Println(i)
	}
}

func (this *Quadrilaterals) printSelf() {
	for _, one := range this.quadrs {
		fmt.Println(one)
	}
}

func (this *Quadrilaterals) findTrisByLine(pointA int, pointB int) (rval []int) {
	tris := this.triangulation.Triangles
	for i := 0; i < len(tris)/3; i++ {
		isContainTheLine := (pointA == tris[i*3+0] && pointB == tris[i*3+1]) ||
			(pointA == tris[i*3+0] && pointB == tris[i*3+2]) ||
			(pointA == tris[i*3+1] && pointB == tris[i*3+0]) ||
			(pointA == tris[i*3+1] && pointB == tris[i*3+2]) ||
			(pointA == tris[i*3+2] && pointB == tris[i*3+0]) ||
			(pointA == tris[i*3+2] && pointB == tris[i*3+1])

		if isContainTheLine {
			rval = append(rval, i)
		}
	}
	return
}

func (this *Quadrilaterals) getQuad(tri1 [3]int, tri2 [3]int) [4]int {
	if tri2 == [3]int{0, 0, 0} {
		return [4]int{tri1[0], tri1[1], tri1[2], tri1[2]}
	}
	tri2ThirdPointIndex := 3
	tri2CommonPointsIndexs := []int{}
	tri1ThirdPointIndex := 3
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			if tri2[i] == tri1[j] {
				tri2ThirdPointIndex -= i
				tri1ThirdPointIndex -= j
				tri2CommonPointsIndexs = append(tri2CommonPointsIndexs, i)
			}
		}
	}
	// fmt.Println(points[tris[tri2*3+tri2ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[0]]], points[tris[tri1*3+tri1ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[1]]])
	return [4]int{tri2[tri2ThirdPointIndex], tri2[tri2CommonPointsIndexs[0]], tri1[tri1ThirdPointIndex], tri2[tri2CommonPointsIndexs[1]]}
}

func (this *Quadrilaterals) isConvexQuad(tri1 int, tri2 int) (rval bool) {
	tris := this.triangulation.Triangles
	points := this.triangulation.Points
	tri2ThirdPointIndex := 3
	tri2CommonPointsIndexs := []int{}
	tri1ThirdPointIndex := 3
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			if tris[tri2*3+i] == tris[tri1*3+j] {
				tri2ThirdPointIndex -= i
				tri1ThirdPointIndex -= j
				tri2CommonPointsIndexs = append(tri2CommonPointsIndexs, i)
			}
		}
	}
	// fmt.Println(points[tris[tri2*3+tri2ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[0]]], points[tris[tri1*3+tri1ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[1]]])
	return isConvexQuad(points[tris[tri2*3+tri2ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[0]]], points[tris[tri1*3+tri1ThirdPointIndex]], points[tris[tri2*3+tri2CommonPointsIndexs[1]]])
}

func Triangulate(Lng, Lat []float64) (*delaunay.Triangulation, error) {
	// Create a set of points
	points := []delaunay.Point{}
	for i := range Lng {
		points = append(points, delaunay.Point{X: Lng[i], Y: Lat[i]})
	}

	// Perform Delaunay triangulation
	triangulation, err := delaunay.Triangulate(points)
	return triangulation, err
}

func triTofile(lng, lat []float64) {
	var sb gobase.BytesBuilder
	sb.AppendStr(`{
  "features": [
    {
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [`)

	for i := 0; i < len(lng)/3; i++ {
		if i > 0 {
			sb.WriteString(",")
		}
		sb.Appendf("[[")
		for j := 0; j < 3; j++ {
			if j > 0 {
				sb.AppendStr(",")
			}
			sb.Appendf("[%.6f,%.4f]", lng[3*i+j], lat[3*i+j])
		}
		sb.Appendf(",[%.6f,%.4f]]]", lng[3*i], lat[3*i])
	}

	sb.WriteString(`] 
      },
      "properties": {
        "fill-opacity": 1,
        "stroke-width": 1
      },
      "type": "Feature"
    }
  ],
  "type": "FeatureCollection"
}`)
	// t.Logf("%s", sb.String())
	sb.SaveToFile("assnet.geojson")
	fmt.Printf("------------\n")
}
func pointTofile(ps []delaunay.Point, filename string) {
	var sb gobase.BytesBuilder
	sb.Appendf("%v\n", "name,L,B")
	for i := 0; i < len(ps); i++ {
		sb.Appendf("%v,%.6f,%.4f \n", i, ps[i].X, ps[i].Y)
	}

	sb.SaveToFile(filename + ".txt")
}

// 计算向量AB和向量BC的叉积
func cross3Product(A, B, C delaunay.Point) float64 {
	AB := delaunay.Point{X: B.X - A.X, Y: B.Y - A.Y}
	BC := delaunay.Point{X: C.X - B.X, Y: C.Y - B.Y}
	// 计算叉积
	return AB.X*BC.Y - AB.Y*BC.X
}

// isConvexQuad检查四边形是否是凸四边形，顶点需要按顺序给出
func isConvexQuad(A, B, C, D delaunay.Point) bool {
	crossAB_BC := cross3Product(A, B, C)
	crossBC_CD := cross3Product(B, C, D)
	crossCD_DA := cross3Product(C, D, A)
	crossDA_AB := cross3Product(D, A, B)
	return (crossAB_BC >= 0 && crossBC_CD >= 0 && crossCD_DA >= 0 && crossDA_AB >= 0) || (crossAB_BC <= 0 && crossBC_CD <= 0 && crossCD_DA <= 0 && crossDA_AB <= 0)
}

// QuadraticLinearInterpolation 实现四边形线性插值
func QuadraticLinearInterpolation(A, B, C, D, innerPoint delaunay.Point) bool {
	crossAB_BC := cross3Product(A, B, C)
	crossBC_CD := cross3Product(B, C, D)
	crossCD_DA := cross3Product(C, D, A)
	crossDA_AB := cross3Product(D, A, B)
	return (crossAB_BC >= 0 && crossBC_CD >= 0 && crossCD_DA >= 0 && crossDA_AB >= 0) || (crossAB_BC <= 0 && crossBC_CD <= 0 && crossCD_DA <= 0 && crossDA_AB <= 0)
}

// LinearInterpolation 实现线性插值
func LinearInterpolation(v0, v1, t float64) float64 {
	// s2.Interpolate(v0, v1, t)
	return (v1-v0)*t + v0
}

// calculateT 计算t
func calculateT(A, B, C, D, P delaunay.Point, AValue, BValue, CValue, DValue float64) float64 {
	// pointTofile([]delaunay.Point{A, B, C, D, P}, "ABCDP")
	x := P.X
	y := P.Y
	DeltaXab := B.X - A.X
	DeltaYab := B.Y - A.Y
	DeltaXdc := C.X - D.X
	DeltaYdc := C.Y - D.Y
	a := DeltaXab*DeltaYdc - DeltaYab*DeltaXdc
	b := x*DeltaYab - x*DeltaYdc + D.Y*DeltaXab + A.X*DeltaYdc + y*DeltaXdc - y*DeltaXab - D.X*DeltaYab - A.Y*DeltaXdc
	c := y*D.X - y*A.X - A.Y*D.X - x*D.Y + A.X*D.Y + A.Y*x
	// fmt.Println(DeltaXab, DeltaYab, DeltaXdc, DeltaYdc, a, b, c)
	Delta := b*b - 4*a*c
	t1 := (-1*b - math.Pow(Delta, 0.5)) / 2 / a
	t2 := (-1*b + math.Pow(Delta, 0.5)) / 2 / a
	var t float64
	if t1 <= 1 && t1 >= 0 {
		t = t1
	} else if t2 <= 1 && t2 >= 0 {
		t = t2
	} else {
		t = -1
	}
	if t != -1 {
		x1 := A.X + t*DeltaXab
		x2 := D.X + t*DeltaXdc
		x1Value := LinearInterpolation(AValue, BValue, t)
		x2Value := LinearInterpolation(DValue, CValue, t)
		tt := (x - x1) / (x2 - x1)
		return LinearInterpolation(x1Value, x2Value, tt)
	} else {
		return calculateTByTri(D, B, C, P, AValue, BValue, CValue)
	}
}

func calculateTByTri(A, B, C, P delaunay.Point, AValue, BValue, CValue float64) float64 {
	// pointTofile([]delaunay.Point{A, B, C, D, P}, "ABCDP")
	x := P.X
	y := P.Y
	DeltaXab := B.X - A.X
	DeltaYab := B.Y - A.Y

	b := x*DeltaYab + C.Y*DeltaXab - y*DeltaXab - C.X*DeltaYab
	c := y*C.X - y*A.X - A.Y*C.X - x*C.Y + A.X*C.Y + A.Y*x
	// fmt.Println(DeltaXab, DeltaYab, DeltaXdc, DeltaYdc, a, b, c)
	t := c / (-1 * b)

	// fmt.Println("t: ", t1, t2)
	x1 := A.X + t*DeltaXab
	// y1 := A.Y + t*DeltaYab
	x2 := C.X
	// y2 := C.Y
	x1Value := LinearInterpolation(AValue, BValue, t)
	x2Value := CValue
	tt := (x - x1) / (x2 - x1)
	return LinearInterpolation(x1Value, x2Value, tt)
}

func crossProduct2(p1, p2 delaunay.Point) float64 {
	return p1.X*p2.Y - p2.X*p1.Y
}

func (this *Quadrilaterals) isPointInTriangle(point delaunay.Point, tri [3]int) bool {
	A := this.triangulation.Points[tri[0]]
	B := this.triangulation.Points[tri[1]]
	C := this.triangulation.Points[tri[2]]
	AB := delaunay.Point{X: B.X - A.X, Y: B.Y - A.Y}
	BC := delaunay.Point{X: C.X - B.X, Y: C.Y - B.Y}
	CA := delaunay.Point{X: A.X - C.X, Y: A.Y - C.Y}

	AP := delaunay.Point{X: point.X - A.X, Y: point.Y - A.Y}
	BP := delaunay.Point{X: point.X - B.X, Y: point.Y - B.Y}
	CP := delaunay.Point{X: point.X - C.X, Y: point.Y - C.Y}

	crossABP := crossProduct2(AB, BP)
	crossBCP := crossProduct2(BC, CP)
	crossCAP := crossProduct2(CA, AP)

	return (crossABP >= 0 && crossBCP >= 0 && crossCAP >= 0) || (crossABP <= 0 && crossBCP <= 0 && crossCAP <= 0)
}
